Clear["Global`*"]
BodePlot[Φ_,f0_,fn_]:=Module[{ω},
GraphicsColumn[{LogLinearPlot[20*Log10[Abs[Φ/.{s->I*ω}]],{ω,10^f0,10^fn},AxesOrigin->{10^f0,0}],
LogLinearPlot[180/Pi*Arg[Φ/.{s->I*ω}],{ω,10^f0,10^fn},AxesOrigin->{10^f0,0}]}]
]
Step[Φ_,tn_]:=Module[{t},Plot[InverseLaplaceTransform[Φ/s,s,t],{t,0,tn}]]

Матрица перехода между связанной и нормальной СК

(Subscript[A, gb])=(({{Cos[ψ], 0, -Sin[ψ]}, {0, 1, 0}, {Sin[ψ], 0, Cos[ψ]}}));

(Subscript[B, gb])=(({{Cos[ϑ], Sin[ϑ], 0}, {-Sin[ϑ], Cos[ϑ], 0}, {0, 0, 1}}));

(Subscript[C, gb])=(({{1, 0, 0}, {0, Cos[γ], Sin[γ]}, {0, -Sin[γ], Cos[γ]}}));

MatrixForm[(Subscript[D, gb])=(Subscript[C, gb]).(Subscript[B, gb]).(Subscript[A, gb])]

MatrixForm[(Subscript[D, bg])=Transpose[(Subscript[D, gb])]]

(Subscript[C, gb])=.;(Subscript[B, gb])=.;(Subscript[A, gb])=.;
(Subscript[D, gb])=ReplaceAll[(Subscript[D, gb]),{ψ->ψ[t],ϑ->ϑ[t],γ->γ[t]}];
(Subscript[D, bg])=ReplaceAll[(Subscript[D, bg]),{ψ->ψ[t],ϑ->ϑ[t],γ->γ[t]}];

Матрица перехода между связанной и скоростной СК

(Subscript[A, ab])=(({{Cos[β], 0, -Sin[β]}, {0, 1, 0}, {Sin[β], 0, Cos[β]}}));

(Subscript[B, ab])=(({{Cos[α], Sin[α], 0}, {-Sin[α], Cos[α], 0}, {0, 0, 1}}));

MatrixForm[(Subscript[D, ab])=(Subscript[B, ab]).(Subscript[A, ab])]

(Subscript[B, ab])=.;(Subscript[A, ab])=.;
(Subscript[D, ab])=ReplaceAll[(Subscript[D, ab]),{α->α[t], β->β[t]}];

Матрица перехода между связанной и полусвязанной СК

(Subscript[D, be])=(({{Cos[α[t]], -Sin[α[t]], 0}, {Sin[α[t]], Cos[α[t]], 0}, {0, 0, 1}}));

Динамические уравнения поступательного движения центра масс ЛА

Уравнения движения центра масс ЛА получим на основании закона изменения количества движения центра масс:

((d(OverVector[Q]))/dt)=(Underscript[Σ, i])(OverVector[(Subscript[F, i])]),
(OverVector[Q])=m(OverVector[V]).

Продифференцируем (OverVector[V]) по правилам дифференцирования вектора, записанного в подвижном базисе:

((d(OverVector[V]))/dt)=((d(OverTilde[V]))/dt)+((OverVector[ω])×(OverVector[V])).

Воздушная скорость:

(Subscript[V, a])={V[t],0,0};

Выразим скорость в связанной СК через воздушную:

MatrixForm[(Subscript[V, b])=(Subscript[D, ab]).(Subscript[V, a])]

Угловая скорость вращения ЛА вокруг центра масс:

(Subscript[ω, b])={(Subscript[ω, x])[t],(Subscript[ω, y])[t],(Subscript[ω, z])[t]};

Сила тяжести в базисе нормальной системы координат:

(Subscript[G, g])={0,-G,0};

Сила тяги двигателя в базисе связанной системы координат:

(Subscript[P, b])={P*Cos[(Subscript[ϕ, p])],P*Sin[(Subscript[ϕ, p])],0};

Равнодействующая аэродинамических сил в базисе скоростной системы:

(Subscript[R, a])={-X,Y,Z};

Запишем уравнение в связанной СК, а затем спроектируем в полусвязанную:

tm=MatrixForm[Simplify[m*(Subscript[D, be]).(D[(Subscript[V, b]),t]+(Subscript[ω, b])\[Cross](Subscript[V, b]))]]==MatrixForm[Expand[Simplify[(Subscript[D, be]).((Subscript[D, gb]).(Subscript[G, g])+(Subscript[D, ab]).(Subscript[R, a])+(Subscript[P, b]))]]/.{P Cos[(Subscript[ϕ, p])] Cos[α[t]]-P Sin[(Subscript[ϕ, p])] Sin[α[t]]->P Cos[(Subscript[ϕ, p])+α[t]],P Cos[α[t]] Sin[(Subscript[ϕ, p])]+P Cos[(Subscript[ϕ, p])] Sin[α[t]]->P Sin[(Subscript[ϕ, p])+α[t]]}];
Column[DYPD=Table[tm[[1,1]][[i]]==tm[[2,1]][[i]],{i,1,3}]]
Clear[tm];

Варьирование уравнения нормальных сил

UKS=DYPD[[2]]

Исходное основное движение - ПГП:
(Subscript[V, 0])=const; (Subscript[ω, x0])=0; (Subscript[ω, y0])=0; (Subscript[ω, z0])=0; (Subscript[β, 0])=0; (Subscript[γ, 0])=0; (Subscript[Z, 0])=0;

IOD=UKS/.{V->(Subscript[V, 0]),β->(Subscript[β, 0]),α->(Subscript[α, 0]),γ->(Subscript[γ, 0]),ϑ->(Subscript[ϑ, 0]),(Subscript[ω, x])->(Subscript[ω, x0]),(Subscript[ω, y])->(Subscript[ω, y0]),P->(Subscript[P, 0]),X->(Subscript[X, 0]),Z->(Subscript[Z, 0])}

IOD=IOD/.{(((Subscript[V, 0])))'[t]->0,(Subscript[β, 0])[t]->0,(Subscript[γ, 0])[t]->0,(Subscript[ω, x0])[t]->0,(Subscript[ω, y0])[t]->0,(Subscript[Z, 0])->0}/.{G Cos[(Subscript[ϑ, 0])[t]] Sin[(Subscript[α, 0])[t]]-G Cos[(Subscript[α, 0])[t]] Sin[(Subscript[ϑ, 0])[t]]->G Cos[(Subscript[α, 0])[t]-(Subscript[ϑ, 0])[t]]}

IOD=IOD/.{(Subscript[α, 0])[t]-(Subscript[ϑ, 0])[t]->(Subscript[θ, 0])[t],
((Subscript[ω, z])[t]-(((Subscript[α, 0])))'[t])->D[(Subscript[θ, 0])[t],t],
-G Cos[(Subscript[α, 0])[t]] Cos[(Subscript[ϑ, 0])[t]]-G Sin[(Subscript[α, 0])[t]] Sin[(Subscript[ϑ, 0])[t]]->-G Cos[(Subscript[α, 0])[t]-(Subscript[ϑ, 0])[t]]}

Исходное возмущенное движение:

IVD=UKS/.{V->(Subscript[V, 1]),β->(Subscript[β, 1]),α->(Subscript[α, 1]),γ->(Subscript[γ, 1]),ϑ->(Subscript[ϑ, 1]),(Subscript[ω, x])->(Subscript[ω, x1]),(Subscript[ω, y])->(Subscript[ω, y1]),P->(Subscript[P, 1]),X->(Subscript[X, 1]),Z->(Subscript[Z, 1])}

IVD=IVD/.{(((Subscript[V, 1])))'[t]->(ΔV)'[t],(Subscript[β, 1])->Δβ,(Subscript[γ, 1])->Δγ}

Пренебрежем вторым порядком малости:

IVD=IVD/.{Cos[Δβ[t]]->1,Sin[Δβ[t]]->0,Cos[Δγ[t]]->1,((Subscript[ω, z])[t]-(((Subscript[α, 1])))'[t])->D[(Subscript[θ, 1])[t],t]} /.
{-G Cos[(Subscript[α, 1])[t]] Cos[(Subscript[ϑ, 1])[t]]-G Sin[(Subscript[α, 1])[t]] Sin[(Subscript[ϑ, 1])[t]]->-G*Cos[(Subscript[α, 1])[t]-(Subscript[ϑ, 1])[t]]}

IVD=IVD/.{G Cos[(Subscript[ϑ, 1])[t]] Sin[(Subscript[α, 1])[t]]-G Cos[(Subscript[α, 1])[t]] Sin[(Subscript[ϑ, 1])[t]]->G Sin[(Subscript[α, 1])[t]-(Subscript[ϑ, 1])[t]]}/.{(Subscript[α, 1])[t]-(Subscript[ϑ, 1])[t]->(Subscript[θ, 1])[t]}

IVD=IVD/.{(Subscript[α, 1])[t]->(Subscript[α, 0])[t]+Δα[t],(Subscript[θ, 1])[t]->(Subscript[θ, 0])[t]+Δθ[t],(Subscript[P, 1])->(Subscript[P, 0])+ΔP,(Subscript[X, 1])->(Subscript[X, 0])+ΔX, ((Δθ)'[t]+(((Subscript[θ, 0])))'[t])->0}

Разложим функции в ряд Тейлора, оставив члены порядка малости не выше первой:

Column[serT={Sin[Δθ[t]+(Subscript[θ, 0])[t]]->Normal[Series[Sin[Δθ[t]+(Subscript[θ, 0])[t]],{Δθ[t],0,1}]],Cos[(Subscript[ϕ, p])+Δα[t]+(Subscript[α, 0])[t]]->Normal[Series[Cos[(Subscript[ϕ, p])+Δα[t]+(Subscript[α, 0])[t]],{Δα[t],0,1}]],
Sin[(Subscript[ϕ, p])+Δα[t]+(Subscript[α, 0])[t]]->Normal[Series[Sin[(Subscript[ϕ, p])+Δα[t]+(Subscript[α, 0])[t]],{Δα[t],0,1}]],
Sin[Δα[t]+(Subscript[α, 0])[t]]->Normal[Series[Sin[Δα[t]+(Subscript[α, 0])[t]],{Δα[t],0,1}]],
Cos[Δα[t]+(Subscript[α, 0])[t]]->Normal[Series[Cos[Δα[t]+(Subscript[α, 0])[t]],{Δα[t],0,1}]],
ΔX->D[X[V],V] ΔV+D[X[α],α] Δα,
ΔP->D[P[V],V] ΔV+D[P[(Subscript[δ, сг])],(Subscript[δ, сг])] (Subscript[Δδ, сг])}];

Column[serT={Cos[Δθ[t]+(Subscript[θ, 0])[t]]->Normal[Series[Cos[Δθ[t]+(Subscript[θ, 0])[t]],{Δθ[t],0,1}]],
Sin[(Subscript[ϕ, p])+Δα[t]+(Subscript[α, 0])[t]]->Normal[Series[Sin[(Subscript[ϕ, p])+Δα[t]+(Subscript[α, 0])[t]],{Δα[t],0,1}]],
ΔP->D[P[V],V] ΔV+D[P[(Subscript[δ, сг])],(Subscript[δ, сг])] (Subscript[Δδ, сг])}]

IVD=IVD/.serT

Возмущенное движение:

VD=IVD[[1]]-IOD[[1]]==FullSimplify[Expand[FullSimplify[IVD[[2]]-IOD[[2]]]]/.{ΔP Δα[t]->0,(Subscript[P, 0]) Δα[t]->0,Δα[t] ΔV->0, Δα[t]ΔV->0,Δα[t] (Subscript[Δδ, сг])->0}] /.{-Cos[(Subscript[θ, 0])[t]]+Cos[(Subscript[α, 0])[t]-(Subscript[ϑ, 0])[t]]->0}

Линеаризованное уравнение

LK=-Δϑ+(Subscript[c, 4]) Δα+(Subscript[e, 2]) ΔV+(Subscript[c, 10]) Δϑ+(Subscript[c, 9]) (Subscript[Δδ, в])+(Subscript[r, 2]) (Subscript[δ, сг])==0;

Анализ устойчивости продольного движения

Критерием допустимости разделения продольного движения на короткопериодическое и длиннопериодическое движения самолета является условие устойчивости короткопериодического движения при (Subscript[с, 2]) > 0.

У первого режима (Subscript[c, 2])=-0.71 ⇒ КПД неустойчиво.
У второго режима (Subscript[c, 2])=118.6 ⇒ КПД устойчиво.

Рассмотрим первый режим

c1=0.208;c2=-0.71;c3=3.72;c4=0.411;c5=0.11;c6=7.48;

Автоматы продольной управляемости дают возможность обеспечить удовлетворительные для ручного управления характеристики устойчивости и управляемости, как динамические, так и статические, для любого, в том числе и неустойчивого в короткопериодическом движении самолета.

Рассчитаем передаточные числа АПУ.

Заданное время регулирования:

Treg=2.5;

Параметр, определяющий нуль передаточной функции (при T0 = 1.5 переходный процесс выходода на заданную перегрузку близок к граничному монотонному):

T0=1.5;

Время регулирования характерное для T0 = 1.5:

Treg0=3.8;

Ω0=(Treg0/Treg)

μ=((3 Ω0-c1-c4-c5)/c3)

CyCgp=((c4 c6)/9.81)

K2=((3 ((Ω0)^2)-c2-c1 c4-c3 c4 μ)/(c3 CyCgp))

MatrixForm[ωz={s+c1+c3 μ,-1,-1}]

MatrixForm[α=Simplify[{CyCgp(c5 s+c2)+CyCgp K2,CyCgp(s+c4),0}]]

MatrixForm[θ={0,0,s}]

MatrixForm[δv={c3,0,0}]

MatrixForm[m1=Transpose[{ωz,α,θ}]]

MatrixForm[m2=Transpose[{ωz,α,δv}]]

W=((Det[m2])/(Det[m1]))

Step[W,10]

BodePlot[W,-2,3]

Рассмотрим второй режим

c1=1.2;c2=118.6;c3=27.55;c4=2.32;c5=0.445;c6=7.28;

Ωa=(Sqrt[c2+c1 c4])

ζa=((c1+c4+c5)/(2 Ωa))

k=5.84 c4 c6

(k/(Ωa))

Т.к.  k/omega>16,  необходим демпфер тангажа.

Рассчитаем параметры демпфера тангажа

μ=((2(Sqrt[c2+c1 c4])-(c1+c4+c5))/c3)

Постоянная времени фильтра высоких частот в цепи сигнала ωz

Tb=(1/c4)

Самолет без демпфера

MatrixForm[ωz={s+c1,-1,-1}]

MatrixForm[α={c5 s+c2,s+c4,0}]

MatrixForm[θ={0,0,s}]

MatrixForm[δv={c3,0,0}]

MatrixForm[m1=Transpose[{ωz,α,θ}]]

MatrixForm[m2=Transpose[{ωz,α,δv}]]

W=((Det[m2])/(Det[m1]))

Step[W,10]

BodePlot[W,-3,3]

Самолет c демпфером

MatrixForm[ωz={s+c1+((c3 μ Tb s)/(Tb s+1)),-1,-1}]

MatrixForm[m1=Transpose[{ωz,α,θ}]]

MatrixForm[m2=Transpose[{ωz,α,δv}]]

W=Simplify[((Det[m2])/(Det[m1]))]

Step[W,10]

BodePlot[W,-3,3]